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Weak Hi physisorption energies present a significant challenge to even the best correlated theoretical many- 
body methods. We use the phaseless auxiliary-field quantum Monte Carlo (AFQMC) method to accurately 
predict the binding energy of Ca + -4H2. Attention has recently focused on this model chemistry to test the 
reliability of electronic structure methods for Hi binding on dispersed alkaline earth metal centers. A modified 
Cholesky decomposition is implemented to realize the Hubbard- Stratonovich transformation efficiently with 
large Gaussian basis sets. We employ the largest correlation-consistent Gaussian type basis sets available, up to 
cc-pCV5Z for Ca, to accurately extrapolate to the complete basis limit. The calculated potential energy curve 
exhibits binding with a double-well structure. 

PACS numbers: 64.70.K-, 71.15.-m, 61.50.Ks, 71.15.Nc. 



I. INTRODUCTION 

Clean energy economy has become an appealing world- 
wide endeavor because of its promise of environmental friend- 
liness and economic security. Developing improved hydro- 
gen fuel storage systems for fuel cell applications and finding 
new technologies for production of hydrogen from renewable 
sources are important components in achieving this goal. A 
major standing challenge is the lack of a method for efficient 
hydrogen storage and retrieval. In developing technologies 
for on-board hydrogen storage systems, the U.S. DOE has set 
a target capacity of 9 wt % by year 20154- Physisorption of hy- 
drogen molecules is one mechanism under consideration for 
use in hydrogen storage. For operation at ambient tempera- 
tures, suitable storage media should have a binding strength of 
~ 0.1 — 0.4 eV per H2^ However, our current understand- 
ing of physisorption processes is still lacking. 2 First-principles 
calculations could potentially play an important role in pro- 
viding insight into the physics of physisorption and identify- 
ing suitable media and mechanisms. Achieving predictive ac- 
curacy for such weakly bound systems has been problematic, 
however, for the usually successful density functional theory 
(DFT) approach and even for explicitly correlated methods. 

Dispersed alkaline-earth metal systems have recently 
shown some promise as H2 storage medial— Attention has 
focused on a simple model chemistry to test the reliability 
of electronic structure methods in predicting the binding of 
H2 on Ca + centers. (This is a simplified model which does 
not take into account the zero-point motion and entropy ef- 
fects, which are important in modeling real hydrogen stor- 
age.) A common thread from these calculations is the pre- 
diction and characterization of a double-well potential energy 
curve (PEC) in the symmetric dissociation of Ca + -4H2. An 
outer van der Waals well is found in some of the correlated 
calculations, but not in the DFT calculations with standard lo- 
cal or semilocal exchange-correlation functionals, 710 An in- 
ner well (Kubas comple x 12 ' 13 ) is found at shorter distances. 
Overbinding is found in local and semilocal DFT calculations, 
while explicitly correlated methods yield conflicting results 
regarding the magnitude of the binding or whether this in- 



ner well is even bound. Disagreement between the many- 
body calculations can be partly attributed to inadequate basis 
set convergence. The earliest calculation, using the second- 
order M0ller-Plesset perturbation method (MP2), predicted an 
inner local minimum that is not bound.- Ohk et al& subse- 
quently used MP2 with larger basis sets with an added correc- 
tion from coupled cluster with singles and doubles and pertur- 
bative triples [CCSD(T)] to find both bound inner and outer 
wells. As we will show in Sec. [Hi] however, the basis sets 
in Ref. [i| were still too small for accurate extrapolation to the 
complete basis set (CBS) limit, due to neglect of (semi)core- 
valence correlation effects. The most recent many-body cal- 
culation used diffusion Monte Carlo (DMC) to find an un- 
bound inner well and a barely bound outer well. 10 The current 
status of theoretical work is clearly unsatisfactory, and the rea- 
sons for this need to be clarified and remedied. 

The primary objectives of our work are to produce an ac- 
curate first-principle PEC of the Ca + -4H2 symmetric dis- 
sociation process and to clarify the key factors that con- 
trol the accuracy for many-body calculations. The phaseless 
auxiliary-field quantum Monte Carlo (AFQMC) methodi 4 ^ 
with standardized Gaussian type orbital (GTO) basis^ here- 
after designated as GTO-AFQMCr^^ is used to compute the 
PEC and to carefully extrapolate the result to the CBS limit. 
AFQMC has been applied to study a wide variety of ma- 
terial systems ^ 18 ' 20 " — Its accuracy has been shown to be 
similar to the gold-standard CCSD(T) method near equilib- 
rium geometries and better for bond-breakingi 17 - 23 i 25 - 26 The 
key features of AFQMC include: (1) low algebraic scaling 
with the system size [O [M 3 — Af 4 ), where M is the number 
of the one-particle basis functions]; (2) ease of implementa- 
tion on massively parallel supercomputers; (3) its wide appli- 
cability as an ab-initio method in condensed matter physics 
and quantum chemistry. In the context of quantum chem- 
istry, GTO-AFQMC uses the identical basis and Hamiltonian 
of other standard methods, such as Hartree-Fock (HF), MP2, 
CCSD(T), and configuration interaction (CI). 

In this paper we implement a significant technical improve- 
ment in GTO-AFQMC, which removes a computational bot- 
tleneck in the preprocessing and initialization step. This allow 
us to use very large basis sets to obtain accurate ground state 
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energies at the CBS limit. This is made possible by the imple- 
mentation of a modified Cholesky decomposition (mCD)22r— 
of the two-body interaction term in the Hamiltonian. Our 
GTO-AFQMC calculation with the aug-cc-pCV5Z basis set 
(827 GTOs) represents the largest many-body calculation in 
this system with GTOs to date. 

The rest of this paper is organized as follows. The imple- 
mentation of mCD together with relevant methodological de- 
tails of AFQMC are presented in Sec.Q]] Accuracy and timing 
illustrations of AFQMC/mCD are also given. GTO-AFQMC 
Ca + -4H2 PEC calculations and extrapolation to the CBS limit 
are presented in Sec. [HI] In Sec.[IV]we discuss possible resid- 
ual sources of error and compare our results to previously pub- 
lished results. We summarize and conclude in Sec.lVl 



n. GTO-AFQMC METHODOLOGICAL IMPROVEMENTS 
AND CALCULATION DETAILS 

After a brief AFQMC overview, methodological improve- 
ments with GTOs are presented. The improvements remove a 
preprocessing bottleneck for computer time and memory, al- 
lowing the use of large GTO basis sets in AFQMC. This is 
achieved in an unbiased manner using mCD. Accuracy and 
timing illustrations are given, and the section concludes with 
the calculation details which are used in Sec. [Ill] 



A. AFQMC ground state projection 

We briefly review the key features of the phaseless AFQMC 
method, which are relevant for the remainder of this section. 
Detailed descriptions of the method can be found in Refs.0- 

M 

AFQMC stochastically evaluates the ground state of a 
many-body Hamiltonian H by means of importance sampled 
random walks in Slater-determinant spaceJ 4 ^ The stochasti- 
cally generated determinants are the samples of the ground- 
state many -body wave function l^o). Although AFQMC is in 
principle an exact method, in practice the sign or complex- 
phase problem arises ; 14 - 30 i 31 which must be controlled us- 
ing some approximate means in order to prevent the expo- 
nential growth of the Monte Carlo variance. This is done 
with the phaseless approximation^ which has been demon- 
strated to yield excellent agreement with both experimen- 
tal and known exact results in a wide variety of molecules 
and extended systems i ' 17 ' 1 8 i 20- £& To date we have developed 
AFQMC methods with two basis sets: (1) planewaves with 
pseudopotentials, which are suitable for studying periodic 
systems j 14 ' 1 8 - 20 - 22 i 24 i 32 (2) GTOs, which are widely employed 
for ab initio quantum chemistry ^I i 2i i -23 1 25 1 26 j n action there 
have been many applications on lattice models {e.g., the Hub- 
bard model) where the phaseless approximation reduces to 
the constrained path approximation 30,33 due to the nature of 
the interaction. These systems typically have correlation en- 
ergy as a significantly higher fraction of their total energy than 
in most molecules and solids. The high accuracy that has 



been achieved in these systems for both energy and correla- 
tion functions is thus encouraging for the prospect of AFQMC 
in real materials. 

The electronic Hamiltonian H = K + V contains one -body 
K = J^ij Kij c l c j kinetic energy and external potential terms 
plus two-body electron-electron interaction terms, 

V = \j2 V ^ klC H c k c i > (1) 

ijkl 

expressed here in a second quantized form. The fermionic cre- 
ation operators c\ are defined on a finite set of M orthonormal 
one-particle basis functions {xi(r)}. 

The AFQMC method projects the ground state wave func- 
tion |^o) from a trial wave function |\& T ), 

e- rA e- T " ■■■e- TA \^> T ) -> |tf„> , (2) 
using a short-time Trotter-Suzuki decomposition 

e -rH = e -rK/2 e -rV e -rK/2 + q ^ (3) 

The input ^> T can be either a single Slater determinant or 
a multi-determinant wave function; the AFQMC projection 
yields a stochastic multi-determinant representation of "fo- 

The one-body projector is straightforward to implement in 
AFQMC: e~ TA '/ 2 acting on a Slater determinant yields an- 
other Slater determinant. To evaluate the two-body projector 
e~ rV , we first decompose V into a sum of squares of one- 
body operators,— 

V = - - "7 + (one-body term) , (4) 

7 

where the minus sign is just a notational convention, and ad- 
ditional one-body terms can arise from reordering the cre- 
ation and destruction operators. We employ the Hubbard- 
Stratonovich (HS) transformation 3 ^ 3 ^ to rewrite the two-body 
projector as a multi-dimensional integral, 

e CV2)rE^ = TT f°° ^l e -^/2 e v^MV . (5) 

This integral over {er 7 } is evaluated stochastically in 
AFQMC. Both Eqs. (f3]) and © are exact in the limit 
t — y 0. We therefore have an exact reformulation of the 
original ground-state projector in terms of one-body projec- 
tion operators {t) 7 } coupled with external auxiliary fields 
{er 7 } which, after integration, recovers the original two-body 
interactions^ 4 - 



B. Modified Cholesky decomposition 

1. A potential computational bottleneck 

The HS decomposition [Eq. (0]i] of the two-body interaction 
term is a preprocessing step, which is done only once at the be- 
ginning of the calculation. The decomposition is not unique, 
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and this flexibility can be exploited to obtain better perfor- 
mance and/or accuracy of the AFQMC calculation. 34 The two- 
body interaction matrix elements Vijui in Eq. (Q~|) are given by 
electron-electron Coulomb repulsion integrals (ERIs) in elec- 
tronic structure calculations: 



Vi 



ijkl 



J dridr-2 X*( r i)xK r 2)' 



ri -r 2 



rXz(ri)Xfe(r2) 



(6) 



In our earlier implementation of GTO-AFQMC, 17 we used the 
straightforward approach of diagonalizing the M 2 x M 2 sym- 
metric ERI supermatrix V^un^u^, where /x = (i, I) and 
v = (j,k) are compound indices. (The dimension M 2 is re- 
duced by a factor of two with GTOs as they are real-valued.) 
With the eigenvalues A 7 and eigenvectors Xj^ » of the real, 
symmetric V^„ supermatrix, the HS one-body operators can 
be written as 



i.l 



(7) 



In practice, only O(M) of the eigenvalues are found to have 
magnitudes greater than ~ 10~ 8 Eh, so the remainder can be 
discarded. The direct diagonalization leads to an 0(M 6 ) scal- 
ing of computer time and 0(M A ) storage. While exact for any 
Hermitian two-body interaction, this approach clearly scales 
poorly with system size. 

The HS decomposition described above is a bottleneck for 
treating large systems. For special choices of the basis, such 
as planewaves, or for model systems, such as the on-site Hub- 
bard model, the two-body interaction is easily written into the 
bilinear HS form of Eq. ©, with O(M) terms. This and 
the diagonalization results above suggest that the information 
content in the two-body term is only O(M) and that more ef- 
ficient HS strategies could be devised for general basis sets in 
electronic structure calculations. 

The underlying problem here is the sheer size of the two- 
body supermatrix, which plagues all ab initio quantum chem- 
istry methods. A number of approaches have been devised to 
reduce the computer time and number of integrals that need 
to be stored. These include density fitting or other auxiliary 
basis methodsr^i resolution of Coulomb operator*^ 3 - and 

We have chosen to implement the mCD to carry out the HS 
decomposition in AFQMC. The method, similar to planewave 
approaches, has a single threshold parameter S which deter- 
mines the maximum error in the Cholesky-expanded ERIs. 
For symmetric, positive semidefinite V^ v supermatrices, this 
error can be reduced to zero, within machine precision, for 
sufficiently small 8. In practice, our calculations use S in the 
range 1CP 6 to 10 .Eh, depending the target statistical accu- 
racy in the AFQMC calculation, resulting in N^d < 7.5M 
Cholesky vectors. The algorithm and its performance are dis- 
cussed next. 



2. Implementation of Cholesky decomposition in AFQMC 

The symmetric, positive semidefinite ERI supermatrix 
is decomposed using a recursive mCD algorithm. 2 ^ 29 Given 
J Cholesky vectors (j = 1...J), can be expressed as 



= V (-J) + a( j ) 



(8) 



where A\fJ is the residual error at the J-th iteration. The 
(J + l)-th Cholesky vector is obtained from 



tJ+1 



A 



(J) 



(9) 



where [v] j indicates the index of the largest diagonal element, 
Al J ) , , , of the residual error matrix in the J-th iteration. 

A key point is that only a single column A^j ; need to be 
computed and stored at any iteration. This iteration is repeated 
until A?, , is less than S. This procedure guarantees that 
all matrix elements of the residual matrix are less than 5: 



VlLV 



V [IV 



A(Ncu) 



< 5 



(10) 



where the number of Cholesky vectors corresponding to this S 
is denoted by Nqb- Damped prescreeningsS, is applied in each 
iteration step to stabilize the mCD factorization. 

Using the -ZVcd Cholesky vectors, the decomposition in 
Eq. becomes 



7=1 \ a 1 \ jk 

^ Ncu 

" O EE i J(ij)^(j,k) C «' 



(11) 



7=1 ijk 

0{5) , 



where the extra one-body operator is also explicitly expressed 
in terms of the Cholesky vectors. The Cholesky vectors trans- 
late directly to the matrix elements of v 1 in Eq. (0): 



(12) 



and the number of auxiliary fields is just given by -/Vcd- 

In this work the Cholesky vectors were generated "on the 
fly" within NWChemii or MPQC.- 4 ^ In GTO-AFQMC, the 
second-quantized expression for the Hamiltonian must be ex- 
pressed with respect to orthogonalized GTOs.— The mCD 
procedure [Eq. ([H])] is carried out in the original GTO basis, 
and the resulting Cholesky vectors are then transformed to the 
orthogonalized basis. 
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In order to produce high-quality Cholesky vectors which 
satisfy the accuracy condition Eq. ( [Tol l, it is imperative that 
the original ERIs used in the mCD recursive procedure be cal- 
culated to very high accuracy and precision. The V^ u super- 
matrix will not be strictly positive semidefinite (within ma- 
chine precision) if there are errors in the calculated ERIs. In 
a test case with M = 180, we observe from direct diago- 
nalization of V^„ that there are some negative eigenvalues 
of O (— 10~ 8 ). In this case, the resulting Cholesky vectors 
would not properly reconstitute all the ERIs to within 8 ac- 
curacy when <5 is driven below 10~ 8 . In another test with 
M ~ 550, errors in the calculated ERIs smaller than lO -8 ^ 
lead to violation of Eq. ([TOj with 8 set to O (lO -6 ) . 



3. GTO-AFQMC/mCD accuracy and timing illustrations 

Table U compares the accuracy of GTO-AFQMC calcula- 
tions using mCD with that using direct diagonalization (DD) 
of the Vjj„ supermatrix. For each 8 shown, we confirmed 
that all V^„ matrix elements are reproduced with error less 
than <5. Within statistical uncertainty, the QMC energies are 
equivalent whether we use DD or mCD with 8 ranging from 
10~ 8 through 10~ 3 For this system, the truncation bias 
from mCD exceeds the targeted statistical error of 2 x 10~ 4 
E h when 8 > 3 x 10~ 3 -E^. Figure Q] plots the error 



TABLE I. GTO-AFQMC total energies for several values of the 
mCD threshold parameter 8 for Ca + -4H 2 (Z = 2.3 A and 
d H _ H = 0.7682 A; see text), using the cc-pVTZ basis (M = 155). 
A fixed Trotter time step r = 0.01-E^ 1 was used for the all the calcu- 
lations. The total energy Edd, obtained from direct diagonalization 
of Vfiu, is presented for comparison; the eigenvalue cutoff is also 
shown. N-y is the corresponding number of auxiliary-fields. A full 
rank, symmetric, positive definite Vy, v matrix would have required 
155 2 = 24025 Cholesky vectors. All energies are in E^. 
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1727 


10~ 6 


-681.43003(17) 


1120 


10" 5 


-681.42988(20) 


850 


HT 4 


-681.42977(18) 


643 


10" 3 


-681.42932(19) 


511 


2 x 10~ 3 


-681.42991(19) 


468 


3 x 10" 3 
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-681.39679(19) 
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5 x 10~ 3 


-681.39609(17) 


379 



-EmCD (8) — -Edd from Table U For comparison, the corre- 
sponding error in the UHF variational energy, computed using 
the same set of Cholesky vectors, is also shown and is seen to 
correlate very well with the AFQMC energies. 
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FIG. 1. (Color online) AFQMC total energy error, E mC D(5) - -Edd, 
as a function of the mCD threshold parameter 5, where Edd is the 
energy obtained from direct diagonalization of the Coulomb matrix 
(energies are given in Table IJl. The width of the gray line indicates 
the statistical uncertainty of Edd- The corresponding error (with re- 
spect to the UHF energy) of the trial wave function variational energy 
is also shown. 




FIG. 2. (Color online) Log-log plots of wall clock times (seconds) 
vs. basis size M, comparing DD (+) and mCD (A) of the ERI su- 
permatrix V^v. The DD is carried out with 8-thread OpenMP (see 
text). The mCD timing is obtained from a locally modified MPQC 45 
program, for fixed 8 = 10 -6 . No multithreading is used in the mCD 
calculations. The dashed (dash-dot) lines are linear regressions. The 
DD slope ~ M 5 ' 7 is consistent with the expected A/ 6 scaling, while 
mCD scales as ~ A/ 3 ' 1 . 



Figure [2] illustrates the relative timing of DD vs. mCD 
procedures. Computations were carried out on 64-bit AMD 
Opteron multi-core processors with speeds ranging from 2.2 
GHz up to 3 GHz.— As expected DD times scale as M 6 . 
The recursive mCD algorithm with prescreening scales only 
as ~ M 3 . For the largest (M = 827) basis reported in this pa- 
per, mCD required less than six hours on a single core, while 
DD would have taken more than 92 days and nearly 1 TB of 
memory on the same computer. 
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Our current mCD implementation still has much room for 
improvement for application to larger systems. The mCD 
for the preprocessing HS decomposition has not been paral- 
lelized. The sparsity of the Cholesky vectors has also not been 
exploited for the actual GTO-AFQMC calculations. The de- 
composition of the necessary V^ v matrix elements scales as 
O (M 3 ) in computer time and the memory required to store 
the Cholesky vectors also currently scales as O (M 3 ). Fully 
exploiting the sparse structure of the Cholesky vectors would 
further reduce the memory requirement to O (M 2 ) asymptot- 
ically for very large molecules^ 



C. Ca + - 4H2 calculation details 

As a model of the Ca + binding site in a hydrogen storage 
system, the PEC for symmetric dissociation of Ca + -4H2 was 
calculated as a function of the Ca + - H2 lateral distance Z (see 
Fig- ID. For each value of Z, the H-H distance c?h-h is opti- 
mized using MP2 with the cc-pCVTZ basis. The UHF wave 
function is used as the trial wave function ^ T . GTO-AFQMC 
calculations were carried out using the correlation-consistent 
core-valence (cc-pCVxZ) basis set family for C a 47 i 48 and cc- 
pVxZ for the H atoms. (This joint basis is subsequently des- 
ignated as "cc-pCVxZ".) For some selected geometries, cal- 
culations were also carried out using a second basis set fam- 
ily denoted "aug-cc-pCVxZ", which comprises the cc-pCVxZ 
functions for Ca and the diffuse aug-cc-pVxZ functions for 
the H atoms. 



a minimum of three basis sets, while the correlation energy 
CBS extrapolation requires at least two. 

We examined errors from the Trotter time step r on basis set 
extrapolation in GTO-AFQMC. With core-valence basis sets, 
the absolute energy can vary significantly. At the 5Z level, for 
example, the error in the absolute energy is as large as 20 mE^ 
for r = 0.01-Eh- Binding energies, however, are less sensitive, 
since the t — > extrapolation slope for a given basis set was 
found to be insensitive to the geometry of the system. Binding 
energies were therefore obtained using a finite time step of 
t = O.OlEr . The validity of this approach was verified by 
computing energy differences obtained with separate r — > 
extrapolations at representative geometries. For the valence- 
only basis sets, the total energies change insignificantly under 
t — y extrapolation, and the calculations are always reported 
atr = O.OLE^ 1 . 

The PEC for the Ca + -4H2 binding energy is given by 



Eb(Z) = E(Z) — -Efrag , 

where Ef lag is the fragment total energy, 

Efia, g = E c& * + 4i?H 2 ■ 



(15) 



(16) 



In our CBS extrapolation it is convenient to further decompose 
the binding energy into its HF and correlation contributions, 



E h (Z)=E^(Z)+E c h °"(Z). 



(17) 



In this work the HF total energies are extrapolated using 
x G {3, 4, 5} in E^ F {Z). The ansatz in Eq. CK implies that 

^corr 
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FIG. 3. An illustration of the Ca + -4H2 model chemistry, containing 
one Ca + surrounded by four hydrogen molecules in a D^h symmetric 
configuration. 

Extensive basis-set extrapolation tests were carried out to 
eliminate the errors from the use of finite GTO basis sets. 
Such tests have not been systematically reported in previous 
studies of the system. The best CBS extrapolation procedure 
separately treats the HF and correlation energies, 



-Ehf(oo) « Euf(x) - Be 

-^corr(oc) ~ E CQrT 



( \ C 

0) - Za ' 



(13) 
(14) 



C'(Z) is a Z-dependent CBS coefficient. Ideally we would 
perform AFQMC calculations at all geometries with two or 
more basis sets to obtain both E^°"'(Z, 00) and C'(Z) di- 
rectly. Such a calculation would be unnecessarily expensive, 
especially at the largest cc-pCV5Z basis level (M = 627). In- 
stead, AFQMC CBS parameters were obtained at several rep- 
resentative geometries (Z) using both the x € {3, 4} and the 
x G {3, 4, 5} series. We then adopt the following strategy to 
parametrize Cafqmc(^) anc ' estimate the final PEC. We com- 
pute the AFQMC finite-basis E^ ori (Z, x) with the cc-pCVTZ 
basis set (x = 3). We also assume that, across all Z values, the 
ratio of correlation energies recovered by AFQMC and MP2 
is approximately constant and given by the parameter p. The 
AFQMC CBS correction term can therefore be approximated 
by scaling the MP2 CBS correction term, 



Ai?™ rr (Z, x, AFQMC) 



pA^ orr (Z,a;,MP2), (19) 



since correlation energy convergence is much slower than 
HF42r— H ere x [ s the correlation consistent basis cardinal 
number. The HF CBS extrapolation requires calculations with 



or, equivalently, 



AFQMC 



(Z)^pCi AP2 (Z) 



(20) 
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The MP2 energies were computed using the cc-pCVTZ and 
cc-pCVQZ basis sets to obtain C' MV2 (Z), which determines 
p and the final CBS estimate for the complete AFQMC PEC. 
We verified that this MP2-assisted approach accurately repro- 
duced direct AFQMC extrapolations at selected geometries to 
within statistical errors. 

Spectroscopic constants associated with the computed PEC 
were obtained from Morse fits, 



E h (Z) = E a 



k 
2^2 



1 - e 



-a(Z-Zo) 



(21) 



where Eq is the well depth minimum, Zq the location of the 
well minimum, and k is the one-dimensional harmonic spring 
constant. 



III. RESULTS 

In modeling the binding of FF molecules onto dispersed 
calcium centers, the convergence with respect to basis set is 
very delicate. This is illustrated by the all-electron Ca + -4H2 
PECs in Fig. @] The solid lines are Morse fits to GTO- 
AFQMC and MP2 binding energies calculated with the cc- 
pCVTZ basis set (see Sec. Ill Cl for computational details). The 
symmetric dissociation PECs are seen to exhibit a double-well 
structure. AFQMC results using larger (cc-pCVQZ and cc- 
pCV5Z) basis sets are shown at a few selected Z. The inner 
and outer wells show dramatically different rates of basis set 
convergence. The inner well does not exhibit binding at the 
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FIG. 4. (Color online) All-electron GTO-AFQMC and MP2 PEC of 
the Ca + -4H2 symmetric dissociation as a function of Z (Ca + - H2 
separation distance) at the cc-pCVTZ basis level. The solid lines are 
from separate Morse fits for the inner and outer regions. Also shown 
at selected Z are GTO-AFQMC results from larger cc-pCVQZ and 
cc-pCV5Z basis sets. Vertical lines at Z = 2.3 A and Z = 4.0 A are 
a guide to the eye. 

cc-pCVTZ level, being ~ 0.2 eV higher than the dissociation 
limit. The outer van der Walls well is bound (~ — 0.1 eV at 
Zq ~ 3.5 A). As the basis size is increased, the inner well 
changes character, becoming ~ OeV with cc-pCVQZ and 
then bound by ~ — 0.1 eV at the cc-pCV5Z level. By con- 
trast, the outer well binding energy is already converged at 
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FIG. 5. (Color online) Basis set convergence of GTO-AFQMC 
Ca + -4H2 total energies for Z = 2.3 A (near the inner well mini- 
mum). Energies are plotted as a function of x~ 3 , where x £ {3, 4, 5} 
is the correlation consistent basis cardinal number. Left panel: 
valence-only cc-pVxZ and aug-cc-pV^Z; right panel: core-valence 
cc-pCVa;Z and aug-cc-pCVxZ. 



cc-pCVQZ basis level. Thus careful extrapolation with the 
basis sets is required in these systems to reach the CBS limit. 

We now show that the basis set must well represent the 
semicore Ca 3s and 3p states for accurate extrapolation. Fig- 
ure [5] plots the total energies of Ca + -4H2 for correlated va- 
lence basis sets (cc-pVi-Z and aug-cc-pVrrZ) and correlated 
core-valence basis sets (cc-pCV^Z and aug-cc-pCVxZ), as a 
function of basis cardinal number [x = (3, 4, 5); see Eq. (fT4bl. 
The calculations are for Z = 2.3 A, which is near the inner 
well minimum. With the core-valence basis sets, the energies 
show linear dependence on x~ 3 , consistent with the ansatz in 
Eq. ( fT4b . while the valence-only series do not. For reference 
and benchmark purposes, selected GTO-AFQMC total ener- 
gies are tabulated in Table HI1 

The importance of proper core-valence treatment is also ev- 
ident in the binding energies, as shown by Fig.|6j where results 
for both the inner- and outer-well regions are shown. For com- 
parison, MP2 results are also shown for the inner-well region. 
The magnitude of core-valence effects is clearly larger in the 
inner-well region. This is due to the Kubas interaction^^i!! 
involving Ca(3d) states. Since the spatial extent of the semi- 
core Ca(3s,3p) is similar to the Ca(3<f) states, core-valence 
effects are magnified in this region. The CBS extrapolation 
lowers the inner well minimum by nearly 0.4 eV compared to 
the cc-pCVTZ basis results. By contrast, the outer well depth 
is lowered by only < 0.1 eV. Even at the cc-pCV5Z level the 
binding energy is still ~ 0.05 eV higher than the CBS limit 
for the inner well, while at the the outer well it has long con- 
verged. Figure|6]also shows that the CBS limit is relatively in- 
sensitive to the use of diffuse "aug" functions for the H atoms. 
These results emphasize that it is necessary to use larger core- 
valence correlation-consistent basis sets in many-body calcu- 
lations for such weakly bound hydrogen storage systems with 
binding-site atoms containing semicore states. 

In Figs.|5]and|6]we see that, with the cc-pCVrrZ basis sets, 
even the total and binding energies follow the x~ 3 scaling to 
a very good degree, which is the form for correlation ener- 
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TABLE II. All electron GTO-AFQMC Ca + -4H2 total energies (in Eh) for two geometries and two correlation-consistent core- valence basis 
sets. "Inner well" corresponds to Z — 2.3 A and g!h-h = 0.7682 A; "outer well" corresponds to Z — 4.0 A and cZh-h = 0.7362 A. The total 
energies of the isolated Ca + and H2 fragments are also shown. The mCD method with 5 — 10 -6 Eh is used, unless otherwise noted. In all 
cases, bias from mCD is much smaller than the AFQMC statistical error, which are on the last two digits and are indicated in parentheses. All 
energies were extrapolated to the r — > limit. CBS indicates the complete basis set extrapolated limit. 



a <5 : 

b 5: 



10" 
10" 



5 E h 







inner well 


outer well 


Ca + 


H 2 


Basis set 


M 


Z = 2.3 A 


Z = 4.0 A 




d H - H = 0.7362 A 



Basis family: cc-pCVa;Z 

x = 3 180 -681.61620(61) -681.62600(71) -676.93411(28) -1.17257(11) 

x = 4 358 -681.83561(74^ -681.84005(76)" -677.13994(45) -1.17384(18) 

x = 5 627 -681.9189(12)" -681.92106(77)" -677.21997(52) -1.17412(15) 

CBS 00 -681.9976(10) -681.99954(84) -677.29517(52) -1.17447(17) 



Basis family: aug-cc-pCVa;Z 

x = 3 252 -681.62437(79) -681.63180(89) 

x = 4 486 -681.8395(11) -681.84151(71) 

x = 5 827 -681.92056(80)^ -681.9213(10) b 

CBS 00 -682.00018(94) -681.99718(99) 



-1.172826(83) 
-1.17397(22) 
-1.17434(12) 
-1.17466(14) 
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FIG. 6. (Color online) Top panels: The binding energy of Ca + -4H2 
system near the inner well minimum (Z — 2.3 A, dn-H = 0.7682 A) 
plotted against (a; -3 ), where x is the correlation consistent basis car- 
dinal number. Bottom panels: The binding energy near the outer well 
minimum (Z = 4.0 A, d H -H = 0.7362 A). 



gies, as shown in Eq. ( [l~4"l i. This is so because the HF en- 
ergy converges rapidly with the basis size, as indicated by 
Eq. ( fl~3l l. For example, the HF energy changes by roughly 
— 4mi?h from cc-pCVTZ to the CBS limit across different Z 
values, in contrast to a change in correlation energy of almost 
— 400m£'h- Similarly, the corresponding HF binding energy 
changes only by about —1 mE^ (—0.03 eV). Thus, using the 
procedure described in Sec. Ill CI the many-body results with 



TABLE III. Binding energies Eh after extrapolation to the CBS 
limit. The two geometries are the same as in Table [TTJ (inner well 
at Z = 2.3 A and outer well at Z = 4.0 A). The contributions to the 
AFQMC binding energy [Eh(Z, 00) on the third row] from HF and 
correlation are shown separately in the first two rows. Energies are 
in eV. AFQMC statistical errors are shown in parentheses. 





Inner well 


Outer well 


E™{Z,oo) 


+0.815 


-0.082 


E C h OII (Z,^) 


-0.954(23) 


-0.072(22) 


E b (Z,oo) 


-0.139(23) 


-0.154(22) 



core-valence correlation-consistent basis sets can be extrapo- 
lated to the CBS limit in a very robust fashion. 

Table [ill] shows the extrapolated binding energy results at 
two representative geometries. HF is qualitatively correct for 
the outer well, capturing more than half of the well depth. For 
the inner well, HF is unbound by a very large amount on the 
scale of interest, by more than 0.8 eV. Thus, consistent with 
the discussion above on the slow convergence of the Kubas 
complex, the binding of the inner well is dominated by elec- 
tron correlation effects, which contribute almost 1 eV to the 
well depth. 

Figure |7] presents the final Ca + -4H2 symmetric PEC from 
GTO-AFQMC calculations, after extrapolation to the CBS 
limit. The GTO-AFQMC binding energies from the aug-cc- 
pCV5Z basis set are shown at Z = 2.3 and 4.0 A to indicate 
the effect of the extrapolation (recall Fig . |4j) . The extrapola- 
tion of the AFQMC results was done with the procedure de- 
scribed in Sec. Ill CI assisted by MP2 CBS corrections. As 
mentioned, this procedure was verified by direct extrapola- 
tions of the AFQMC results at two representative geometries 
using the x £ {3,4,5} series, as shown in Table Hill which 
gave completely consistent results. In addition, at multiple 
other geometries, extrapolations of the AFQMC results using 
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• MP2 CBS 

X AFQMC (CBS) 

T AFQMC / aug-cc-pCV5Z 




Z (angstrom) 

FIG. 7. (Color online) The PEC of the Ca + -4H2 symmetric dissoci- 
ation extrapolated to the CBS limit. The PEC is shown as a function 
of Ca + - H2 separation distance, Z. Morse fits to the AFQMC points 
are shown as solid curves. The gray shading provides an estimate 
of the PEC uncertainties. Dotted lines shows the continuation of the 
outer well PEC into the small Z region. Morse fits to the MP2 points 
are shown as dashed curves. The extrapolation scheme is discussed 
Sec. Ill CI For comparison, binding energies computed using very 
large aug-cc-pCV5Z basis sets are also shown at Z = 2.3 and 4.0 A. 



TABLE IV. Spectroscopic constants associated with the AFQMC in- 
ner and outer wells depicted in Fig. [7] These parameters were ob- 
tained by fitting a Morse curve to the two wells separately. The 
error bars shown in parentheses below include both the fitting and 
AFQMC statistical uncertainties. 



Quantity 


Units 


Inner well 


Outer well 


Well depth minimum Eq 


eV 


-0.178(16) 


-0.1566(71) 


Equilibrium distance Zq 


A 


2.205(17) 


3.375(41) 


Spring constant k 


eV/A 2 


13.0(22) 


0.342(37) 



the x G {3, 4} series were done as further confirmation, which 
lead to statistically indistinguishable results at the CBS limit. 

The corresponding MP2 results at the CBS limit are also 
shown. MP2 is seen to represent the double-well structure 
reasonably well (Figs. |4]and|7). However, there are some sig- 
nificant deviations from AFQMC. The MP2 inner well is more 
shallow: at the largest basis set (cc-pCV5Z), the MP2 inner- 
well binding energy is ~ 0.05 eV smaller than AFQMC. The 
MP2 outer well, by contrast, is deeper and broader than that 
of AFQMC. This is perhaps not surprising, given the much 
stronger effect of electron correlation in the inner well region. 
We comment more on the MP2 results in Sec. lIVI 

The solid lines with gray shading in Fig.UJrepresent Morse 
fits of the GTO- AFQMC binding energies at the CBS limit, 
where the width of the shading indicates the uncertainties in 
the fits. As mentioned, the two regions have separate Morse 
fits. The fitting uncertainty in the inner well is < 0.02 eV, 
and < 0.01 eV for the outer well. Spectroscopic constants 
corresponding to the Morse fitted GTO- AFQMC CBS results 
are shown in TablellVI The two well minima have comparable 
well depths within the AFQMC statistical accuracy. 



IV. DISCUSSION 

In the previous section, the GTO-AFQMC Ca + -4H 2 PEC 
was shown to have a double-well structure with weak binding 
of four H2 molecules. After extrapolation to the CBS limit, 
binding energies of —0.178(16) and —0.1566(71) eV were 
found at the respective minima of the inner and outer wells, 
essentially the same within statistical errors. The crossover 
between the two wells occurs at Z ~ 2.5 A. Below we first 
comment further on possible sources of uncertainty. 

A possible source of error in the PEC is the basis set super- 
position error (BSSE). For a finite basis, when two or more 
fragments are brought together, the resulting "molecule" en- 
joys additional degrees of freedom not present in the isolated 
fragments. This BSSE results in the artificial enhancement of 
the predicted binding energy. Counterpoise (CP) corrections 53 
represent an attempt to reduce the BSSE. In the Ca + -4H2 sys- 
tem, the effect of the CP correction on HF energies is negli- 
gible even at the relatively small cc-pCVTZ basis level. For 
many-body calculations, we investigated the effect of the CP 
correction within MP2. The binding energy at the CBS limit 
was changed by ~ 0.02 eV for the inner well, while the cor- 
rection in the outer well is negligible. This estimate indicates 
that the effect of CP correction is within the statistical error of 
GTO-AFQMC and does not change our conclusions. 

A second and related possible source of error is the ba- 
sis set convergence, as discussed extensively in the previous 
two sections. This is a system which is particularly demand- 
ing in terms of reaching the CBS limit, as we have already 
shown. The consistency of our various cross-checks suggest 
that the errors from basis set extrapolation are captured in 
the indicated uncertainties in Fig. [7] and in Table [TV] After 
we completed this work we were made aware of more recent 
core-valence basis sets by Iron and co-workers^. They have 
been employed to determine spectroscopic constants for quan- 
tum defect calculations in calcium. 55 ' 56 We have tested these 
new basis sets, with and without the additional d functions, 
using MP2 on the Ca + -4H 2 binding energy at Z = 2.3 A. 
We found that the binding energy results were consistent with 
those from the cc-pCVxZ basis sets which we have employed 
in the present paper. For example, the binding energy at the 
5Z level agreed to within 0.003 eV. 

A third possible source of error is that our AFQMC method 
is not exact. The phaseless approximation is made to control 
the sign/phase problem, and as a result there is a systematic 
bias. The ground state energies calculated from the method 
is not guaranteed to be a variational upper bound. As men- 
tioned, for a variety of benchmarks and applications, the sys- 
tematic bias is shown to be very small, consistently reaching 
the level of accuracy of CCSD(T). For bond-stretching and 
bond-breaking, the method is shown to be more accurate than 
CCSD(T)] 17 i 23 i 25 - 26 For the present systems, internal checks 
by varying the form of the trial wave function indicate that the 
results are very robust. 

We next compare our results with previous calculations. As 
discussed in the previous section, MP2 is seen to describe the 
system quite well. Due to its nonperturbative nature, AFQMC 
recovers more correlation energy compared to MP2. For ex- 
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ample, at the CBS limit, the AFQMC Ca + -4H 2 total ener- 
gies are ~ 60 mE^ lower at all geometries. Cancellation of 
errors greatly reduces the discrepancies in the corresponding 
binding energies. The high symmetry and absence of near- 
degeneracy in the Ca + -4H2 system also makes it easier for 
MP2 to perform well. The fourfold symmetric Ca + -4H2 sys- 
tem is in a half-filled, closed-shell configuration. Open-shell 
configurations and near degeneracies in other cases would be 
more challenging. 

The CCSD(T)/CBS results of Ohk and co-workers^ also 
predict binding for the both the inner and outer wells. How- 
ever, Ohk et al. used valence-only correlation consistent basis 
sets, and only up to cc-pVQZ to find the CBS limit. As shown 
in the previous section, extrapolations with valence-only ba- 
sis sets is problematic. Moreover, their crossover position is 
at Z ~ 2.3 A, while our AFQMC and MP2 crossovers are at 
Z-2.5 A. 

In comparing with the DMC results of Bajdich et al.^ we 
note that in general DMC is a highly accurate many-body 
method. Their calculation predicts an inner well that is not 
bound (Eq ~ 0.06 eV) and an outer well that is weakly bound 
(Eq ~ — 0.07eV). Their crossover position is more consistent 
with our PEC than with the results of Ohk et They used 
Z = 4.6 A as representative of the unbound fragments. Our 
result suggests that at Z = 4.6 A the Ca + -4H2 is still weakly 

bound with E h 0.08 eV. Correcting the DMC PEC by 

this amount increases the DMC outer well binding strength to 
agree with AFQMC. The DMC inner well remains essentially 
unbound. The DMC calculations fixed g?h-h = 0.77 A, but 
our test calculations with AFQMC varying c?h-h show that 
this only has a small effect on the binding energy, < 0.01 eV. 
Thus, the most likely cause for the discrepancy on the inner 
well would appear to be the fixed-node error in DMC. The use 
of only triple-zeta quality basis sets in the trial wave functions 
of the DMC calculations may have contributed to increase this 
error. 



V. SUMMARY 

The phaseless auxiliary-field quantum Monte Carlo method 
was used to accurately predict the binding energy of 
Ca + -4H2, a model chemistry that has recently been used 
to test the reliability of electronic structure methods for 



H2 binding on dispersed alkaline earth centers. A modi- 
fied Cholesky decomposition is implemented to realize the 
Hubbard-Stratonovich transformation efficiently in AFQMC 
with large basis sets, which removes a memory and compu- 
tational bottleneck. We employ the largest correlation con- 
sistent Gaussian-type basis sets available, up to cc-pCV5Z 
for Ca, to accurately extrapolate to the complete basis 
limit. The resulting potential energy curve exhibits bind- 
ing with a double-well structure, with nearly equal binding 
energy minima of ~ — 0.18eV. We showed that an accu- 
rate description of the inner well Kubas complex requires 
the use of the correlation-consistent core-valence basis set 
series cc-pCVxZ. While the model's binding energy of 
~ —0.04 eV/H2 falls short, in itself, of the targeted optimum 
design value by more than a factor of three, the results are en- 
couraging for further study of larger, more realistic models as 
potential candidates for hydrogen storage media. 

The results in this paper demonstrate that the phaseless 
AFQMC method can accurately treat the weakly bound sys- 
tems of interest for hydrogen storage. This is consistent with 
earlier results on a variety of materials systems. As GTO- 
AFQMC continues to improve, it can be expected to treat 
larger, more realistic, nanoscale size systems, with the largest 
available GTO basis sets. We hope that the latest results will 
encourage the integration of the GTO-AFQMC method with 
quantum chemistry approaches, as a component of the com- 
munity's efforts for improving energy technologies. 
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